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SYNTHESIZING COHERENT CORRELATION SUMS AT ONE OR 
MULTIPLE CARRIER FREQUENCIES USING CORRELATION SUMS 
CALCULATED AT A COARSE SET OP FREQUENCIES 
FIELD OF THE INVENTION 

5 The present invention relates to signal processing and, more particularly, to techniques 
for synthesizing coherent correlation sums at a fine frequency resolution using coherent cor- 
relation sums that are calculated at a coarse frequency resolution. 

BACKGROUND OF THE INVENTION 

m In general, when a carrier-modulated signal is received by a receiver, the carrier frequency 
BIO of the signal is shifted by an unknown amount, due to clock drift at the receiver and/or a 
W Doppler shift due to motion. Further, the known signal may be delayed by an unknown 
amount of time in the course of traveling from the source to the receiver. The signal that is 
ftj received at the receiver is generally referred to as the received signal, and it typically includes 
O noise. For many applications, one needs to determine the shift of the carrier frequency to a 
15 high degree of accuracy in order to be able to estimate some of the signal parameters such 
as the signal delay or the data that are being communicated through the signal. 

In one approach for estimating the carrier frequency, known as coherent processing, 
multiple candidate carrier frequencies are examined. For each candidate frequency, In Phase 
("I") and Quadrature ("Q") correlation sums are calculated, followed by an evaluation of 
20 the sum of squares P + Q^. The candidate frequency that results in the largest sum of 
squares is selected as the estimated carrier frequency. If the delay is an unknown quantity, 
then a two-dimensional search is performed to simultaneously estimate the delay and carrier 
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frequency. 

A significant disadvantage to the above approach is that the computational expense in 
performing the calculations for the various candidate frequencies can be prohibitively high. 
This is because the number of candidate frequencies that have to be examined must increase 
5 in proportion to the duration of the signal and therefore in proportion to the number N of 
available samples of the received signal. Consequently, the overall computational effort is of 
the order of at least N'^. 

In another approach, the signal is broken into multiple short blocks, each block is pro- 
5; cessed coherently as above, and then the values of P+Q^ calculated from the different blocks 
qIO are added. This approach, known as non-coherent processing, can work with a smaller num- 
ry ber of candidate frequencies, compared to coherent processing. On the other hand, it is much 
more sensitive to noise. Therefore^ in order to maintain a constant level of performance, more 
m signal samples are required and the overall computational effort remains high. 
□ Based on the foregoing, there is a clear need for a Technique for synthesizing coherent 

15 correlation sums at a fine frequency resolution using coherent correlation sums that are 
calculated at a coarse frequency resolution that requires less computational expense. 

PROBLEM STATEMENT 
The general problem of estimating an unknown shift in the carrier frequency of a received 
signal is as follows. Given a sequence of N real numbers, denoted by xi, . . . , ^rjv, the following 
20 expression is to be evaluated (exactly or approximately) : 

C{f)^j2xke'^^f''^, (1) 
fc=i 
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for a plurality of frequencies / in a given range 

Herein, j is used to denote the square root of —1, that is, j = Also, A is a known 

constant that represents the spacing in time between consecutive numbers Xf^. Note that 
5 C{f) is complex- valued. The real and imaginary parts of C{f) are known as In Phase (I) 
and Quadrature (Q) correlation sums, respectively. 

In a variant of the problem, it is desired to maximize |C'(/)p, commonly referred to as 
S: the "ambiguity function," over the given frequency range. This maximization problem is 
Co typically solved by evaluating over a densely-spaced grid of candidate values of /, 

fillO and selecting that value of / for which \C{f) \^ is largest. This procedure finds the maximum 
L the chosen grid. However, when the grid spacing is small and the function C{f) is 

ni smooth, the results of this search provide a very close approximation to the search over the 
C3 entire range [Fi, F2] of frequencies. In another variant, one may seek to identify values of / 
for which exceeds a certain threshold. 

15 In a further variant of the problem, the numbers Xk depend, in a known manner, on 
an additional unknown parameter a, and are denoted by Xfc(a) to make this dependence 
explicit. In this variant, the problem consists of calculating the correlation sum 

C{f,a)^f:xkia)e'^^^'''', (2) 

k=l 

for various values of the pair (/, a) or of maximizing C{f^a) over a given range of (/, <j) 
20 pairs. 
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In another variant, the parameter a can take integer values within a given range M — 
{0, 1, 2, . . . , M - 1} and Xk[a) is of the form 

for some known sequences yk and Si, resulting in 

, C(/,a)^X:?/^^^-^'"''^'^- (3) 

^ k=i 

In a special case of the problem one may wish to calculate the correlation sum for a 

S single choice or for very few choices of /, but do this in a computationally efficient manner 

0 that makes use of correlation sums that have already been computed for frequencies that 

^ diflFer from the frequencies that are of interest. This situation can arise in a number of 

""10 ways. In some cases, the signal is broken into blocks and correlation sums are automatically 

f\ computed at the receiver, e.g,, by a bank of analog or digital correlators, at a preset set 

^1 of frequencies. One may want to make computationally efficient use of the available block 

correlation sums to synthesize the desired correlation sums over the entire duration of the 

signal. In some cases, the raw data may be unavailable, so that such a synthesis of long 

15 correlation sums on the basis of short correlation sums is the only option available. In other 

cases, an algorithm may have already computed correlation sums for certain values of /, and 

of particular interest are the correlation sums for a nearby value of /. 

As another variant, one may want to compute a correlation sum of a more general form, 

such as 

20 C{f,o) = f:y,5,(a)e^-^•^^^+^•^^ (4) 
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or 

C{f,a)^j:y,s,.y-^^'^-'^'\ (5) 
where 9k is a known phase term that changes slowly with time. 

ORIGINS OF THE PROBLEM 

The above stated problem arises during the acquisition of a carrier-modulated signal 
that has been distorted by an unknown shift of its carrier frequency. In this subsection, a 
summary description of this particular context is provided. 

A receiver receives a signal of the form 

y(t) = a5(t;a)cos(27r(/o + r)t + eW + 0*)+^W, te[0,T]. 

Here, t stands for time, T is the total duration of the signal, /o is a known carrier frequency, /* 
is an unknown shift of the carrier frequency, e{t) is an additional known phase term, (j)* is the 
typically unknown phase of the modulating signal, and w{t) represents noise. Furthermore, 
a is a usually unknown constant which is proportional to the square root of the power of 
the received signal, and s{t; a) stands for the baseband form of the received signal, in the 
absence of noise. Finally, a is an unknown parameter that determines the exact form of the 
signal in the absence of noise. The form of the dependence on a is assumed to be known. 
That is, if a were known, then the value of s{t; a) would be known for every t. 

In a digital communications context, s{t; a) can be a transmitted waveform, and a can 
represent the data bits that are being communicated through this waveform. In wireless 
and satellite communications, including Global Positioning Systems (GPS), a can represent 
an unknown time delay from the transmitter to the receiver. By letting a be the ratio of 
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the delay to the intersample time A, and if this ratio is integer, one obtains that s{t; a) is 
actually of the form s{t — a* A). 

The variable /* represents an unknown frequency shift, which is known to belong to 
certain range [^1,^2]- This frequency shift can arise in various ways. A representative 
5 but not exclusive list includes receiver motion, transmitter motion, or a clock skew at the 
transmitter or receiver. A frequency shift due to motion is known as a Doppler shift. 

The additional but known phase term d{t) can arise in a situation where the transmitter 
or receiver's velocity is changing in a predictable manner. For a specific context, consider 
^ a GPS receiver that is operating on data that have been collected over a 1 second interval. 
SlO During that interval, a satellite that transmits a GPS signal moves considerably, and the 
fy satellite motion manifests itself as a Doppler shift of the carrier frequency. In addition, 
because the satellite does not move on a straight line relative to a stationary receiver but 
rather follows an elliptical orbit, the velocity of the satellite can change appreciably during 
□ the 1 second interval. Such a change in velocity (acceleration) can be modeled as a change 
15 in the carrier frequency or equivalently as a slowly changing additional phase term 6{t), 
Because the satellite orbit is known, the additional phase term d{t) is also known. 

In some contexts, the signal is immediately mixed at the receiver from the carrier fre- 
quency to an intermediate frequency, in which case the notation /o in the above formula can 
be taken to represent the intermediate frequency. In some contexts, the signal is immediately 
20 mixed directly to baseband, in which case /o can be taken equal to zero. 

In some cases, the received signal y{t) is analog and is subsequently filtered and sampled. 
In other cases, the received signal is obtained from an intermediary device in filtered and 
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sampled form. In either case, the result is a discrete sequence of samples of the form 

Vk - Sk{a) cos(27r(/o + /*)fcA + + (^*) + wj,. 

Here, y^y ^^(a), 9^, and Wk reflect the value of the functions y{t)j s{t;a)^ 9{t) and w{t)^ or of 
filtered versions of these functions, at time A;A, where A is the length of time that elapses 
5 between two consecutive samples. Without loss of generality, it is assumed here that the 
time at which the first sample is taken is A, and that T is an integer multiple of A. Let N 
be the total number of samples, which is given by 

T 

A coherent receiver operates as follows. It calculates In-Phase and Quadrature correlation 
10 sums of the form 

N 

Kf^ <^)^Yj ykSk{(y) cos(27r/A;A + 9k), 
k=i 

and 

N 

for various pairs (/, a), and searches for a pair that results in the largest value of the ambi- 
15 guity function 

One now defines = ykSk{<y), and C{f,a) = I{f,cr) + jQ(/, a). One then observes that 
coherent identification of the carrier frequency for a particular candidate value of a coincides 
with the problem of maximizing a)\^, with C{f, a) defined in Eq. (4), that is, 
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When a is the delay in the reception of the signal, that is, iisit] a) is of the form 5(t-aA), 
then $k{<y) is of the form Skicj) — Sk^^. It follows that coherent identification of the carrier 
frequency for a particular candidate value of a coincides with the problem of maximizing 
|C(/,a)p, with C{f,a) defined in Eq, (5), that is. 

When acquiring a signal, one typically needs to perform a two-dimensional search over 
a set of candidate pairs (/, cr). If the ambiguity function jC(/,<7)p is large for a particular 
(/, cr), this is taken as evidence that this particular pair provides a close approximation to the 
true frequency shift and parameter of the received signal, and the signal is said to be detected 
or acquired. It is common to combine this approach with a thresholding scheme whereby 
acquisition is deemed to be successful if the ambiguity function exceeds some threshold. 

In some cases, the received signal is analog but is not sampled. The processing of the 
signal can be carried out in analog form, using for instance analog correlators implemented in 
hardware. In such cases, all of our previous discussion remains applicable, provided that the 
sums are replaced throughout by integrals. For ease of exposition, our discussion will focus 
on the case where the signal is sampled and sums are calculated. However, one skilled in the 
art would have no difficulty applying this invention to the case where sums are replaced by 
integrals. 

TRADITIONAL METHODS 

In this section, the methods that are traditionally used for carrier frequency estimation 
are discussed, and their shortcomings are indicated. Before doing so, the number of candidate 
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frequencies that need to be examined are discussed. 

THE NUMBER OF CANDIDATE FREQUENCIES 
Since there are infinitely many different frequencies in the range [Fi, F2], this range has to 
be discretized so that C(/, a) is only computed for a finite number of candidate frequencies 
5 /, namely those that belong to a suitably constructed grid of frequencies. 

Let us take the spacing between successive grid points to be a number A/. One then 
computes C(/,cr) only for those / of the form f = Fi + bAf, where b is an integer taking 
O values in the range & 0, 1, . . . , 5. The number B is chosen so that Ft + BAj > F2 and is 
S approximately equal to 

With this choice, there is a guarantee that at least one grid point will be within A//2 from 
m the true frequency /*. 

rU When a grid frequency / is used in place of the true frequency /*, the phase term e^""^^^^ 

^ will be different from e^""-?-^*^^. The difference must be kept small enough in order to maintain 
15 coherence. For this, it is necessary and sufficient to choose the grid-spacing A/ so that A/T 
is smaller than some user-selected threshold By setting 

A/T - ~, (6) 
n 



it is obtained that the number B of frequency grid points is approximately 

B ^ ^1^3. = niF2 - Fi)T = n{F2 - Fi)AN. 
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For any given choice of n, it is important to note that the number of frequency grid points 
increases in proportion to the number N of data samples. 

A reasonable choice is to let n ^ 4, which results in a maximum phase error (loss of 
coherence) smaller than 7r/4, although other choices are also appropriate. When using a 
larger value for n, accuracy is somewhat improved, but computational requirements increase 
because of the larger number B of frequency grid points. 

THE STANDARD METHOD 
The most common approach is to evaluate the expression for C(/, a) separately for each 
one of the B frequency grid points. If one counts one multiplication and one addition of 
complex numbers as one operation, this approach requires about NB operations. Since B 
increases in proportion to iV, the computational requirements of this approach are of the 
order of A^^. 

If is to be evaluated for M different values of a, the above procedure has to be 

repeated M times, for a total computational effort of MBN. 

To put this in perspective, consider the context where 1 second of a received GPS signal, 
sampled at 4MHz, is searched. With a frequency shift uncertainty of IKHz, and with n ^ 4, 
one has 5^ 4, 000, and = 4 x 10^. If one searches for M ^ 1000 possible values of the 
signal delay, the computational requirements amount to 16 x 10^^, just to acquire the signal 
from a single satellite, which is prohibitively large. 

In the further special case where the data Xf. are of the form Xk = ykSk~a^ some savings 
are possible. One recognizes that C(/, a) is the convolution of two sequences, and that the 
convolution only needs to be computed for M different choices of a. One can divide the 
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sequences into blocks of length M, calculate convolutions block-by-block and add them at 
the end. The convolution of two sequences of length M can be calculated by forming the 
transform of the two sequences, using a Fast Fourier Transform (FFT) algorithm, multiplying 
the transforms, and then taking an inverse transform. This takes ScMlogM arithmetic 
5 operations, where cM log M is the running time of a single FFT or inverse FFT. 

The constant c depends on the implementation of the algorithm, but can be assumed to 
be less than 10. The overall complexity, summing over N/M blocks, and over B frequencies 
is ScNBlogM. Here and throughout the rest of this document, all logarithms will be taken 
with respect to base 2. 

LO In all cases, and because B increases linearly with N, the overall computational effort 
grows as 

A TRANSFORM-BASED METHOD 
There is an alternative to the standard method, which is applicable to the case where the 
additional time-dependent phase term 9k is absent, and which is the subject of US patent 
15 4,701,934, entitled, "Method of Doppler Searching in A Digital GPS Receiver" by inventor 
Steven C. Jasper. One notes, that for any fixed a, the expression for C{f,a) 

A;=l 

is the same as the definition of the discrete Fourier transform of the sequence Xk{(T). Thus, 
one can take the FFT of the sequence Xk{cr) and obtain C{f,a) for all frequencies that are 
20 integer multiples of 1/T. In this case, and with the choice A/ = 'i-/{nT), as in Eq. (6), 
one actually needs the transform at a finer frequency resolution, finer by a factor of n. (For 
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examplej if n — 4, the frequency resolution needs to be finer by a factor of 4.) This resolution 
can be obtained by taking the original sequence, of length AT, and appending (n- 1)7V zeroes 
at the end, to obtain a sequence of length niV, and then computing the transform of the 
longer sequence. The overall complexity of this method is that of forming the FFT of a 
5 sequence of length nN^ which is cnN\og{nN), 

For the case of multiple (say, M) candidate signal parameters <7, the process has to be 
repeated M times, resulting in complexity of cnMN\og{nN). 

In the special case where Xk{(y) = ykSk~<j, no further simplification appears possible and 
fn the complexity remains cnMNlog(nN)j which can also be prohibitively large. Furthermore, 
yflO after performing this method^ one does not have a computationally convenient way to adjust 
Ij the sum to another nearby frequency. 

5 NONCOHERENT PROCESSING 

Another common way of avoiding the need for a large number of candidate frequencies 
is to process the data record non-coherently. A typical approach is the following. Break the 
15 data record Xi, . , . , into a number N/Nq (assumed integer) of blocks, consisting of A^o data 
points each. For each data block, use one of the preceding methods to obtain an ambiguity 
function \Ci{f,a)\^ based on the data in the ith block. One then combines noncoherently 
the results from the various blocks, to obtain an overall ambiguity function given by 

N/No 

20 This method can be advantageous because coherence only needs to be maintained over 
intervals of length A^o-^ (the length of a block), as opposed to iVA. Hence, the number 
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of candidate frequencies, and the corresponding computational complexity is reduced by a 
factor of N/Nq, In the absence of noise, this method delivers results comparable to the 
coherent method with much less computation. 

The disadvantage of this method lies in the fact that it is much more sensitive to noise. 
5 In many settings, coherent processing is optimal, as far as its statistical properties are con- 
cerned. Non-coherent processing suffers a so-called "squaring loss" that results in reduced 
sensitivity. This loss becomes more pronounced when the signal is weak and a long data 
record is needed, 

m The invention described in this application aims at reducing computational requirements 

CClO for the case of long data records. Long data records are needed mainly when the signal 
is weak and the signal-to-noise ratio (SNR) is low. In this context, the noise sensitivity 
pi of non-coherent processing is very pronounced. For the non-coherent method to achieve a 

fy sensitivity comparable to that of the coherent one, the data set would have to be much 

ly 

^ larger, and there would be no net decrease in computational requirements. 

15 SUMMARY OF THE INVENTION 

Techniques are provided for synthesizing a long coherent I and Q correlation sum or in- 
tegral at a particular frequency by synthetically combining a sequence of shorter correlation 
sums or integrals at the same or different frequency. Techniques are also provided for acquir- 
ing a carrier-modulated signal with an unknown shift of the carrier frequency, and possibly 
20 some additional unknown signal parameters. These techniques involve synthesizing coher- 
ent correlation sums or integrals at a fine frequency resolution, using coherent correlation 
sums or integrals that are calculated at a coarse frequency resolution. This approach allows 
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for coherent processing of the received signal over an arbitrarily long time interval, while 
avoiding the excessive computational requirements of traditional methods. 
The method provides an approximation D{f,a) to the correlation sum 

k=l 

5 for various candidate values of the pair (/, a), and uses these approximate values to acquire 
the signal, i.e., identify the true value of the carrier frequency and of signal parameter a. 
In the special case where a is a delay parameter, and Sk{o) is of the form Sk-c, further 
O efficiency improvements are obtained by using Fast Fourier Transform techniques. A similar 

2 approximation is provided in the case where it is desired to approximate correlation integrals 
ftllO instead of correlation sums. 

H The central idea is to synthesize long coherent correlation sums, at a fine set of frequencies, 

3 by combining short coherent correlation sums at a coarser set of frequencies. Specifically, 
% if we split the index set {1, . . . , A^} into blocks {/cj, . . . , h+i - 1} and if we calculate short 

correlation sums 

for some frequency /', and for each block {h, h+i - 1}, then we can synthesize a very 
good approximation D(/,(t) of the long correlation sum C{f,a) at a nearby frequency /, by 
letting 

i 

20 The resulting approximation Z)(/, a) is of the form 
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where 4)f{k) is a time-dependent phase such that \2i:(l)f{k)A - 27r//cA - 6^] is small for 
all k. This results in D{f,a) being a close approximation of C{f,a). The time-dependent 
phase 4>f{k) has a particular piecewise linear form, and is linear inside each block. In some 
embodiments^ the method includes a procedure that selects certain parameters and a specific 
5 choice of ^f[k) that optimizes the overall running time. 

PERFORMANCE LOSSES 

With proper choice of the algorithm parameters and of the time-dependent phase <t>f{k), 
iO this method does not suffer from coherence loss any more than the classical approaches do. 
5 A justification of this statement follows. 

{Ao Recall that in classical methods for frequency search, one chooses a frequency grid spacing 
, " A/ so that A/ — 1/nT — l/(nA^A), where n is a user-selected parameter. This ensures the 
CP following property. For every possible frequency /*, there is a frequency / on the selected 
grid, for which C{f^a) is calculated, and for which \ f - /*| < l/(2nA^A), resulting in 

\27rfkA - 27r/*A;A| < — , for all k, 

15 This method will guarantee the following analogous property. For every possible frequency 
/*, there is a frequency / on the selected grid, for which D{f, a) is calculated, and for which 

\27r(j)f{k)A - 2^rkA - Ok\ < — , for all L 

Thus, if one chooses n in the method to be twice as large as the value of n used in a classical 
method, the method provides an approximation of the true phase 27rpkA + that has the 
20 same guarantees as those provided by classical methods. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

The present invention is illustrated by way of example, and not by way of limitation, 
in the figures of the accompanying drawings and in which like reference numerals refer to 
similar elements and in which: 

FIG. lA is a block diagram that illustrates a system overview for processing a signal; 
FIG. IB is a flowchart that illustrates the operational overview of a technique for syn- 
thesizing coherent sums and for using them for the purpose of carrier estimation and signal 
acquisition. 

FIG. 2 A and FIG. 2B are flowcharts that provide an alternative iUustration of a technique 
lo for estimating the carrier frequency of a received signal; 

FIG. 2C is a flowchart that illustrates the details of the process for choosing a particular 
instantiation of the method and its associated parameters; 
FIG. 3 illustrates a structure of a coarse-grained grid; 

FIG. 4 is a flowchart depicting the details of the process for calculating the block-level 
15 correlations; 

FIG. 5 is a flowchart depicting the details of the process for generating the reference 
signal; 

FIG. 6 is a flowchart depicting the details of the process for combining the block-level 
correlations to obtain the correlation function D{f,a); 
20 FIG. 7 illustrates the frequencies of the form + fg comprising the set JF; and 

FIG. 8 is a block diagram illustrating a computer system on which embodiments of the 
invention may be implemented. 
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DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 

Techniques are provided for synthesizing coherent correlation sums or integrals at a fine 
frequency resolution, using coherent correlation sums or integrals that are calculated at a 
coarse frequency resolution. This approach allows for coherent processing of the received 
5 signal over an arbitrarily long time interval, while avoiding the excessive computational 
requirements of traditional methods. In the following description, for the purposes of expla- 
nation, numerous specific details are set forth in order to provide a thorough understanding 
of the present invention. It will be apparent, however, to one skilled in the art that the 
present invention may be practiced without these specific details. In other instances, well- 
10 known structures and devices are shown in block diagram form in order to avoid unnecessarily 
obscuring the present invention. 

FUNCTIONAL OVERVIEW 
FIG. lA is a block diagram that illustrates a system overview for processing a signal. 
The system comprises a signal source S, a receiver H, and a server D. By way of example, the 
15 signal source is a machine that produces an analog signal. The analog signal is known, except 
that it depends on an unknown signal parameter a in a known manner. The analog signal is 
transmitted to receiver H. The signal that is received at H is herein referred to as a "received 
signal" . In some embodiments the unknown signal parameter a is an unknown delay in time 
from the time the analog signal leaves signal source S and is received at receiver H. In one 
20 embodiment of the invention, H converts the analog signal into a sequence of numbers by 
filtering and then sampling the received signal. The sampled received signal is herein referred 
to as "sampled data" or simply "data." In one embodiment of the invention, H transmits 
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the sampled data to server D for processing. Further, for the purpose of explanation, assume 
that the known signal's carrier frequency is modulated by an unknown frequency shift, which 
may, for example, be due to a clock drift at the receiver, and /or relative motion of the receiver 
with respect to the signal source. Thus, it is assumed that the sampled data received at 
5 server D has an unknown carrier frequency due to the unknown frequency shift, and that 
the frequency shift belongs to a known frequency range [^1,^2]- 

The present invention is equally applicable to the case where correlation integrals (instead 
of correlation sums) are to be computed, possibly in hardware, using an analog received 
signal, without first sampling the received signal. For ease of exposition, we describe an 
tlO embodiment for the case where calculations are based on sampled data. However, one skilled 
in the art would have no difficulty in making the necessary adjustments for the case of an 
analog received signal. 

In FIG. IB as explained herein, the method receives certain input parameters 100, which 
consist of the following. N is the length of the data record (number of samples) to be 
15 processed. A is the real time elapsed between consecutive samples in the data record, Fi 
and F2 specify the frequency range to be searched; in particular, the set of frequencies / of 
interest is specified by Fi < / < F2. The parameter n is a positive integer that specifies a 
tolerance parameter for phase errors. Given n, the algorithm will ensure that phase errors 
due to approximations do not exceed 27r/n radians. In particular, the method will examine 
20 frequencies that belong to a grid, where the distance between consecutive grid points is 
approximately l/{nNA) Hz. M is a positive integer that specifies the number of different 
candidate signal parameters a that will be examined. DELAY is a binary parameter. If 
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DELAY is TRUE, the candidate parameters a are integer and the reference signals Sk{c^) are 
of the form s^-a, so that a represents time delay. If DELAY is FALSE, then no assumption 
on the structure of the Sk{o') can be made. In either case, and without loss of generality, one 
can index the candidate signal parameters a so that they range from 1 to M. The notation 
5 — {1, . . . , M} is used to refer to the set of candidate a. Finally, c is a positive number 
that reflects the computational complexity of performing Fast Fourier Transforms (FFTs). In 
particular, the complexity of performing an FFT or an inverse FFT on a (possibly complex) 
sequence of length n will be estimated to be cn log n. 

FIG. IB is a flowchart that illustrates the operational overview of a technique for synthe- 
^0 sizing coherent sums. The process of block 200 is described in greater detail in FIG. 2C. At 
fy block 200 of FIG. IB, a particular instantiation of the method is selected, so as to optimize 
^ its computational requirements. Block 200 receives input parameters 100. The instantiation 
!f ! of the method is specified by certain algorithm parameters 250. These algorithm parameters 
Ps include the length TVq blocks into which the data are broken, the number L of such 

15 blocks, the number H of different coarse-grained frequencies that will be considered while 
processing each block, and the set of these frequencies /i, . . . , /i?. The latter frequencies form 
a grid with constant spacing between consecutive grid points, The last algorithm parameters 
are two binary variables, BLOCK-FFT and COMB-FFT. If BLOCK-FFT is TRUE, then 
an FFT-based algorithm is to be used during the process at block 400 of processing each 
20 block. The details of process at block 400 are described in greater detail in FIG. 4 herein. If 
COMB-FFT is TRUE, then an FFT-based algorithm is to be used by the process at block 
800 of combining the results from the various blocks. The details of process at block 800 is 
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described in further detail in FIG. 6 herein. 

The data yi, . . . , yjv to be processed, shown in block 300, are fed to a process at block 400 
that calculates for each block for / — 1, . , . , for each frequency fh in the set /i, . . . , /^^^^ 
and for each tr in the set {1, . . . , M}, the block-level correlations at block 700, according 
5 to the formula 

INo 

To carry out this task, the process at block 400 submits requests to a process at block 600 
that generates the required values 500 of the reference signals sj^{a). The process at block 
600 is described in greater detail in FIG. 5 herein. 

10 The block-level correlations are complex-valued. Their real and imaginary parts are 
known as I and Q correlation sums. Accordingly, a block-level correlation will also be 
referred to as a "pair of I and Q correlation values" or as a "pair of / and Q correlation 
sums" or as a "pair of I and Q correlation integrals." 

The process at block 800 first constructs a finite set !F of equispaced fine-grained fre- 

15 quencies to be examined, of cardinality at most LH. These frequencies are constructed by 
refining the frequency grid /i, . . . , /if of the algorithm parameters 250 by a factor of L, that 
is, by reducing the spacing between consecutive grid points by a factor of L, while omit- 
ting some of the frequencies that fall outside the frequency range [Fl,F2]. A more detailed 
description of the construction of is provided later. 

20 The process at block 800 takes the block-level correlations ui^h{<^) as inputs and com- 
bines them, by forming weighted sums of the correlation sums ui^h{o'), to produce for each 
frequency / in the set and for each candidate signal parameter (T in ^ {1, . . . , M}, an 
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approximation 850 D{f,a) to the correlation sum C{f,a), so that 

D{f,a) « C{f,a) = ^xM<^)e'^^^'^^^''- 

k=l 

These approximate correlation sums are herein also referred to as weighted sums of 
correlation functions and are squared by the process at block 880 to obtain the approximate 
ambiguity function 900 \D{f, c7)p, for every / G and every a ^ M. Once this approximate 
ambiguity function becomes available for every / € ^ and every a e M, a final process at 
block 950 is used to select a pair (/, a) at which the approximate ambiguity is maximized. 

In one embodiment, the maximizing pair (/, a) is selected when a single best estimate 
of an unknown frequency and signal parameter is to be reported. However, different em- 
bodiments are possible. In some embodiments, the process may quit after process at block 
800 and feed the approximate correlation sums D{f, a) to some other method for further 
processing. In other embodiments, the process may quit after process at block 880 and 
feed the approximate ambiguity function \D{f,a)\'^ to some other method for further pro- 
cessing. Such other methods may employ additional side information, or may perform a 
finer-grained search by interpolating the function D{f, a) as described in U.S. Patent Ap- 
plication No. , entitled " Extracting Fine-Tuned Estimates from Correlation 

Functions Evaluated at a Limited Number of Values " by Anant Sahai, John Tsitsiklis, Ste- 
fano Casadei, Andrew Chou, Benjamin Van Roy and Jesse Robert Stone, filed on the same 
day herewith, Attorney Docket No. 60021-0013. Finally, in some embodiments, the process 
at block 950 may select all pairs (/, a) that exceed a certain threshold. Or the process at 
block 950 may select a fixed number of pairs (/, a) with the highest ambiguity values, for 

Attorney Docket: 60021-0012 -21- 



further processing. 

FIG. 2 A and FIG. 2B are block diagrams that provide an alternative illustration of a 
technique for estimating the carrier frequency of a received signal. At block 202 of FIG. 2Aj 
a frequency range of interest is determined. A frequency range of interest is the range of 
5 frequencies within which it is assumed that the true carrier frequency lies. In one embodiment 
of the invention, the frequency range of interest represents a pre-selected frequency interval 
around the carrier frequency of the known analog signal produced by signal source S. 

At block 204, the frequency range is divided into multiple sets of frequency intervals, 
generally known as frequency bins. In one embodiment, two sets of frequency bins are 
'^"^0 used. The first set of frequency bins comprises a set of coarse grain frequency bins while 
the second set of frequency bins comprises a set of fine grain frequency bins. According to 
certain embodiments of the invention, the number of coarse grain frequency bins is preset to 
a constant which is independent of the length of the sampled data. In other embodiments, 
the number of coarse grain frequency bins is chosen to be proportional to the square root 
15 of the length of the sampled data. In addition, in certain embodiments, the selection of the 
number of coarse grain frequency bins may be based upon a user-selected signal-to-noise 
ratio tolerance. 

At block 206, the sampled data is divided into a set of data blocks. In one embodiment, 
the set of data blocks corresponds to the set of coarse grain frequency bins, in the sense 
20 that the length of the data blocks is chosen by taking into account the size of the coarse 
grain frequency bins. At block 208, a range of hypothesized signal parameters such as delay 
values are determined by estimating the range of possible delay in time that is measured 
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from the time the signal is transmitted from the signal source S until the time that the signal 
is received at receiver H. At block 210, for each coarse grain frequency bin, a trial value of 
frequency is selected. In some embodiments the trial value is the midpoint of a frequency 
bin. The trial value of frequency is herein referred to as a "trial At block 212, control is 
5 passed to block 220 of FIG. 2B. 

At block 220 of FIG. 2B, for each data block, using the selected trial /, the In Phase 
and Quadrature correlation sums ( "/ and Q correlation sums" ) are calculated over the range 
of hypothesized delay values. At block 222, for each fine grain frequency bin and for each 
2 hypothesized delay value, the / and Q correlation sums are suitably weighted and then 
ffiO summed over all the data blocks, without re-calculating the / and Q correlation sums, 
fy At block 224, if I and Q represent the summation of the I and Q sums respectively, then 

1^ for each fine gain frequency bin and for each hypothesized delay value, the magnitude of 
the vector (J, Q) is calculated. At block 226, the fine grain frequency bin that contains the 
o highest magnitude value is determined. The frequency value corresponding to the frequency 
15 bin that has the highest magnitude value is estimated to be the carrier frequency of the 
received signal, which includes the Doppler shift, if any. 

SELECTION OF AN INSTANTIATION OF THE METHOD 
FIG. 2C is a flowchart that illustrates the details of the process at 200 of FIG. IB for 
choosing a particular instantiation of the method and its associated parameters. The process 
20 begins at block 228 of FIG. 20 by generating a set of candidate block lengths iVo. In one 
embodiment, the set of candidate block lengths consists of all positive integers that are 
powers of 2 and that are smaller than N, In alternative embodiments the set of candidate 
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iVo could be any set consisting of integers in the range {2, . . . , iV - 1}. 
At block 229 for each candidate block length Nq^ calculate 

H = \n{F2 - Fi)NoA] , 

which will be the number of frequencies that will be considered for individual blocks and 

which will be the corresponding number of blocks. 

The next steps will rely on the following estimates of the computational complexity 
5 for different instantiations of the processes at block 400 and block 800 of FIG. IB, These 
H estimates are justified later, when the processes at block 400 and at block 800 are described 
IJlO in more detail herein in FIG. 4 and FIG. 6 respectively. 

If block-level correlations are to be calculated by process at block 400 in a default manner, 
y the computational effort of process at block 400 of FIG. IB is estimated as 

t MHN, 

If block-level correlations are to be calculated by process at block 400 of FIG. IB using 
15 an FFT-based method, the computational effort of process at block 400 is estimated as 

3c^H{No + M)logiNo-^M) 

iVo 

This option is available only if Skic) is of the form Sk-a, that is, if the input DELAY is 
TRUE. 

If the combining of the block-level correlations is to be carried out by process at block 
20 800 of FIG. IB in a default manner, the computational effort of process at block 800 is 
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estimated as 

If the combining of the block- level correlations is to be carried out by process at block 800 
using an FFT-based method, the computational effort of process at block 800 is estimated 
5 as 

At block 230 of FIG. 2C, it is determined whether DELAY is FALSE or TRUE. 
If DELAY^FALSE, the process at block 240 of FIG. 2C is used, which does the following. 
For each candidate value of A^^o and associated value of H, it evaluates 



„,^A^, /nN\ 



a ,,,,,, . rMHN^ rr..^ . f»^\ 



and selects that candidate value of A^o for which the above expression is smallest. The above 
1 expression consists of the estimated complexity of calculating the block-level correlations 
^ in the default manner plus the estimated complexity of the better of the two alternative 

available methods of combining the block-level correlations. Furthermore, if the selected 
15 value of Nq, and the associated value of if, satisfy ^^^^ < cniJM-— log then the 

binary variable COMB-FFT is set to FALSE. Otherwise, COMB-FFT is set to TRUE. 
If DELAY=TRUE, the process at block 250 is used, which does the following. For each 

candidate value of iVo and associated value of H, it evaluates 



mm 



[mHN, 3c£i/(JVo + M)log(Aro + M)|+min|^^^^ cnffM^log^^jj 



20 and selects that candidate value of TVq for which the above expression is smallest. The above 
expression consists of the estimated complexity of calculating the block-level correlations 
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using the better of the two alternative available methods plus the estimated complexity 
of combining the block-level correlations using the better of the two alternative available 
methods. Furthermore, if the selected value of A^o? ^nd the associated value of H, satisfy 
MHN < 3c^H{No + M)\og{No + M), then the binary variable BLOCK-FFT is set to 
5 FALSE. Otherwise, BLOCK-FFT is set to TRUE. Finally, if the selected value of No, and 
the associated value of satisfy ^^^^§t- < cnHM^ log j then the binary variable 
COMB-FFT is set to FALSE. Otherwise, COMB-FFT is set to TRUE. 

At the conclusion of process at block 240 or at block 250, a particular value for Nq and 
an associated value of H and L have been selected. Then, the process at block 260 chooses 



ffflO the coarse-grained grid of frequencies /ij . . . , /j/ by letting 

J Because of the way that H and fn are chosen, the highest frequency on the grid is at 

m least — l/(2A^An). This is because 

H>n{F2-Fi)NoA 

15 and 



2iVoAn ' 'No An - ' 2NoAn ' N(,An' 

Since the distance between consecutive grid points is I/INqAu), it follows that for every 

/ G [Fi, F2] there is some fh which is within l/(2iVoAn) from /. 

The structure of this grid is illustrated in FIG. 3. In that figure, and for the purpose of 
20 illustration, it is assumed that Fi = 5000, F2 = 9200, and NqAu = 1/1000. This leads to 
iJ = 5 and /i = 5500, /a = 6500, fs = 7500, U = 8500, h = 9500. 
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CALCULATION OF BLOCK-LEVEL CORRELATIONS 
FIG. 4 is a flowchart depicting the details of process at block 400 of FIG. IB for calculating 
the block-level correlations. For each block / G {1, . . . , L}, for each coarse-grained frequency 
/a, G {1, . . . , F}j and for each a e — {1, . . . , M}, this process calculates the block-level 
5 correlation sums 

INo 

k={l~l)No+l 

The process first examines at block 410 of FIG. 4 the algorithm parameter DELAY. If the 
algorithm parameter DELAY is FALSE, the process receives at block 420 s^a), . . . , S]^{a) 
^ from the reference signal generator, and calculates at block 430 the desired Ui^h{(^) by straight- 
OlO forward multiplications and summations. The total computational effort of this step is esti- 
'^'^ mated as 

^ If the algorithm parameter DELAY is TRUE, the process receives at block 440 S-m+i , . . . , 

I from the reference signal generator and calculates for each block / € {1, . . , , L}, for each fre- 
15 quency fh,he{l,..., H}, and for each a € = {1, . . . , M}, 

INq 

after choosing at block 450 between one of two alternative methods. Thus, at block 450, it 
is determined whether BLOCK-FFT is TRUE or FALSE. 

If BLOCK-FFT is FALSE, it proceeds using at block 460 straightforward multiplications 
20 and summations, and the total computational effort of this step is estimated as 

MHN. 

Attorney Docket: 60021-0012 -27- 



Alternatively, if BLOCK-PFT is TRUE, it uses a process at block 470 that relies on the 
fact that for any fixed I and h, the quantities ui^h{o') correspond to a convolution of the 
sequence s^i^i^No-m+u • - • , sino, of length N^ + M, with the sequence yf^e'^''3h{k-{i~-i)No)A^ 
k = {l~l)No+lj INo, of length iVo. This convolution is calculated by padding both sequences 
5 with zeroes at the end, so that their length becomes the same and equal to Nq + M, taking 
the FFT of the two sequences, multiplying the FFTs, and then taking the inverse FFT. 

The complexity of computing a single convolution is estimated as 3c{No+M) \og{No+M)^ 
where the factor of 3 is due to the fact that there is a total of 3 FFTs or inverse FFTs. There 
2 are LH convolutions to be computed (one for each (/, h) pair) and the total computational 
ailO effort of this step is estimated as 

H 3c^H{No + M) log(iVo + M), 

Alternative embodiments of the calculation at block 400 of FIG. IB for block- level cor- 
^ relations are possible, especially when the reference signal has a special structure. For 
example, when processing GPS signals, the reference signal is a periodic PRN (pseudoran- 
15 dom noise) code, modulated by a slow binary signal. This structure is exploited in U.S. 

Patent Application No, , entitled "SIGNAL ACQUISITION USING DATA 

BIT INFORMATION" by Anant Sahai, Wallace Mann, Andrew Chou and Benjamin Van 
Roy, filed on the same day herewith, Attorney Docket No. 60021-0011. 

Alternative embodiments are also possible in the case where the received signal is not 
20 sampled, but is processed in analog form. In that case, the block-level correlations ui^h{a) 
can be calculated using one or more analog correlators, implemented in hardware. 
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GENERATING THE REFERENCE SIGNAL 
FIG. 5 is a flowchart depicting the details of the process at block 600 of FIG. IB for gen- 
erating the reference signal. At block 510 of FIG. 5, it is determined whether the parameter 
DELAY is TRUE or FALSE. If algorithm parameter DELAY is FALSE, then at block 520 of 
5 FIG. 5 the reference signal generator stores in random access memory a sequence of values 
Si{a)^ . . . ^ s^{a)j for every a E M and supplies these values as needed by process at block 
400 of FIG. IB. If algorithm parameter DELAY is TRUE, then at block 530 the reference 
signal generator stores in its memory a sequence of values 5_m+ij • • • , and supplies these 
^0 values as needed by process at block 400. In particular, for every k and a, Skicr) ~ Sk-a is set. 
In alternative embodiments, the reference signal is not read from memory but is generated 

? 5 ; 

^ on the fly by a computational algorithm. 

. In alternative embodiments, the reference signal is calculated by combining information 

Cfl that is stored in local memory with other information that is externally provided. For 

m 

1^ instance, if the reference signal is a GPS signal, it can be formed by modulating a locally 
15 stored PRN code by a slowly varying binary signal consisting of '^navigation bits." These 
navigation bits are not a priori known, but can be obtained in real time from a commercial 
provider or from a differential GPS station via a communication network, such as the internet. 

COMBINING THE BLOCK-CORRELATIONS 
FIG. 6 is a flowchart depicting the details of the process at 800 of FIG. IB for combining 
20 the block-level correlations to obtain the correlation function cr), for every a e M and 
for every / in a certain set T of fine-grained frequencies that is constructed in the beginning 
of process at block 800. 
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At block 810 of FIG. 6 a set of frequencies /i, . . . , /l, is constructed using the formula 

1 1 q-l 1 ^ , 

2niVoA ^ 2LnA^oA L nNoA' ^~ ' 

Note that these frequencies are spaced l/(nLA^oA) apart, which is a factor of L smaller than 
the spacing l/(niVoA) of the coarse-grained frequencies /i, * . . , /if. Subsequently, at block 
5 820 the process forms the set of all frequencies of the form 

fh + fq, 

i where h ranges over the set (1, . . . , F} and q ranges over the set {1, . . . , L}. Then, at block 
1 830, the above set is pruned by eliminating those frequencies + fq that are larger than 
J F2 + (I/uNqAL), Let be the resulting set, which is the set of fine-grained frequencies to 
^10 be considered. Its cardinality is denoted by B and satisfies B < LH, 
I The smallest frequency in that set is 



2A^oAn 2A^oAn 2LiVoAn 2LiVoAn* 
The largest frequency of the form fh + fq is 

fH+fr = Fi + + (H-1)— + +^~~^ ^ 

JH ^ JL ^ 1 ^ o AT". A ^ ~ ^ AT AT . A ~ O T AT A ^ ^ 



2NoAn ' ^A^oAn 2nNoA 2LNoAn L nNoA 

1^ " 2nNoAL' 

Because of the choice of the spacing, which is l/(nA^oAL), it follows that for every /* E 
[Fx, F2]y there exists some fh + fq ^ ^ which is within l/(2nNoAL) of 

The frequencies of the form fh + fg comprising the set ^ are illustrated in FIG. 7, In 
FIG. 7, and for the purpose of illustration, it is assumed that Fx ^ 5000, F2 = 9200, 
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No An — 1/1000, and L — L In particular, H = 5 and the spacing between consecutive 
fine-grained frequencies is l/(nLiVoA) = 250. The fine-grained frequencies are indicated by 
the heavy dots on the frequency axis. For reference, the coarse-grained frequencies fh are 
also indicated. The lowest fine-grained frequency is Fi + l/{nLNoA) = 5125 and the highest 
5 is 9375. The cardinality B of is 18, which is somewhat smaller than LH — 20. 
The rest of this process at block 800 calculates 

Cl for every fh + fq ^ The calculation uses as input 840 the numbers ui^h{o) (one such 
number for each / G {1, . . . , L}, each G {1, . . . , H}, and each a G M) which are the output 

fHio of the process at block 400 of FIG. IB. At Block 845, it is determined whether the binary 
algorithm parameter COMB-FFT has been set by process at block 200 of FIG. IB to TRUE 

1 or FALSE. 

% If COMB-FFT is FALSE, the process at block 850 that is used to calculate D{fh + fq, a] 

just implements the above displayed formula. This is done for M values of a, and B values of 
15 fh + fq^ For each {fh + fq, o-) pair, there are L ^ N/Nq complex additions and multiplications 
to be performed. Thus, the complexity of the process at block 850 can be estimated as 

No 

Typically, one has B/N ^ H/Nq, or S « NH/Nq, and the above complexity estimate 
becomes 

20 
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If COMB-FFT is TRUE, an alternative process at block 860 is used. It is observed that 
D(fh + fq,o-) can be rewritten as 

1=1 

Viewed as a function of /gNoA, this is the discrete Fourier transform of the sequence 
5 tx,^(fT)e2''J'('-^)^''^°^+^'^('-i^^o, / = l,...,L. Thus, one can take the FFT of this sequence 
and obtain D{fh + /g, o) for all fq such that /^A^o A is an integer multiple of the fundamental 
frequency which is 1/L. That is the same as computing D{fh + iq,o) for all /g that are in- 
j teger multiples of l/(LiVoA). However, the frequencies /, of interest are spaced 1/ {nLNoA) 
§ apart, and therefore, the frequency resolution has to be refined by a factor of n. 
IJlO The above is accomplished by the process at block 860 which does the following. Take 
\1 the original sequence ui,hicj)e'^''^^^-^^^''^°^^''v-im , / = 1, . . . , L, of length L, append {n - 1)L 
O zeroes at the end, to obtain a sequence of length nL, and then compute the FFT of the 
51 longer sequence. For each cr € M and for each fh, the transform at the required resolution 

; VI 

U is computed with cnL log(nL)L operations. Since L » N/Nq, the complexity of this process 
15 is estimated as 

AT , /nN\ 
cn/fM- log (-). 

APPROXIMATION GUARANTEES 
In this subsection, it is explained why the method described in the above has the desired 
approximation guarantees. For simphcity, let us consider the case where the additional phase 
20 term 9k is absent. The quantity approximated by the method is given by 
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The approximation constructed is of the form 

For k in the ^th block, one has 

Mk) = fh{k'{l- l)No) + (/ - + f,)No = (A + f,)k - IS -{I- l))iVo. 

5 Indeed, the first term fh{k - (/ - l)A^o) is introduced when computing the block-level corre- 
lations 

INo 

' k={l-l)No+l 

The second term (l — l)(fh + fqjNo is introduced by the combining formula 

10 For every possible frequency /* in the range [Fi, F2], and because the frequencies fh + fq 
are spaced l/{nNQLA) apart, there is some f = fh + fq such that 

If - ri < — - — ^ 

' - 2nNoLA 

Further, it follows that 

Ifk-Mk)] < \r-fh-Uk + \f,\{k-{l-l)No) 

1 

~~ nA' 

It follows that for every possible frequency /*, there is a fine-grained frequency f = fh + fq, 
for which D{f, a) is calculated, and for which 

|27r<^/(fc)A - 2Trf*kA\ < — , for all k. 
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When the time- varying phase 9k is present, there is an additional phase error equal to the 
maximum change of 6k in the course of a block, When 6k changes slowly, and provided that 
the block length is not too large, the error due to this additional phase term is negligible. 

ALTERNATIVE EMBODIMENTS 
The embodiments as described herein are easily modified to deal with certain changes 
in the objective of the computation or to accommodate certain preferences and other con- 
siderations. A number of alternative embodiments of the various steps are discussed in the 
sequeL 

SELECTION OF COARSE AND FINE-GRAINED FREQUENCY GRIDS 

There are many ways of selecting the grid of coarse frequencies and the grid T of 
fine-grained frequencies. In some embodiments, the density of the grid need not be L 
times more than the density of the initial grid /i, . . . In some embodiments, one or 
both of these grids can be chosen to have non-uniform spacing between consecutive grid 
points. In some embodiments, the frequency grid points are selected ahead of time by the 
user, whereas in other embodiments the frequency grids are selected by the algorithm taking 
into account various input parameters. Finally, in some embodiments, one may choose the 
offsets between the coarse and the fine grids so that the coarse grid is a subset of the fine 
grid. 

CHOOSING THE FREQUENCY GRID SPACINGS 

In one embodiment, the spacing of the coarse and of the fine-grained frequency grids is 
indirectly specified by the input parameter n that reflects a tolerance for phase approximation 
errors. In alternative embodiments, the various grid spacings are directly specified by the 
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user. In other embodimentSj one uses a formula that quantifies performance losses as a 
function of the spacing of the various grids and selects the grid spacing by evaluating the 
tradeoff between performance losses and computational complexity. 
PARTITIONING THE FREQUENCY RANGE [^1,^2]. 

In some embodiments, one can start by partitioning the frequency range [fi,F2] into 
disjoint or overlapping frequency subintervals, and apply the method separately for each one 
of these subintervals. In some instances, the subintervals can be chosen small enough so that 
H — \ and a single coarse-grained frequency will suffice for each subinterval. That single 
frequency can be chosen to be the center of the subinterval, or can be chosen differently. 
SELECTION OF ALGORITHM PARAMETERS 

The selection of a particular method during process at block 200 of FIG. IB, can be 
made using alternative complexity criteria. For example, if certain FFTs (of the reference 
signal or of the data sequence) are available at no extra computation cost, (e.g., if they are 
calculated by some other part of a broader algorithm), then the corresponding computational 
cost should not be included in the cost of process at block 400. Furthermore, if different 
embodiments of the processes at block 400 and at block 800 of FIG. IB are employed, 
then their computational costs have to be estimated and used in selecting the most efficient 
method. 

The particular complexity criteria used in the embodiment of process at block 200 of FIG. 
IB can be replaced by other estimates of the run times of the various blocks, depending on 
the desired degree of precision. For example, the cardinality of the set was estimated as 
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LH. An alternative embodiment estimates the cardinality of T as 

\n{F2 - Fi)iVAl , 

which can be more accurate when F2 - Fi is small. 

One embodiment chooses between FFT-based and default methods, and sets the binary 
variables BLOCK-FFT and COMB-FFT based on computational complexity estimates. In 
alternative embodiments, these choices are preset by the user. 

In alternative embodiments, the method is applied in contexts where the binary input 
variable DELAY is always TRUE, or always FALSE. In such cases, this variable need not 
be treated as an input of the algorithm, and the method can be accordingly simplified by 
omitting those alternatives that do not apply. 

One embodiment provides a specific method for choosing the block size Aq. In alternative 
embodiments the block size can be preset by the user. Or the block size can be selected 
according to another set of criteria. In some cases, there may be a restriction on the minimum 
allowed block size. In some cases, there may be a restriction on the maximum allowed block 
size. In some cases, the block length may be restricted to be approximately equal to the 
square root of A^. 
MISSING DATA 

Certain embodiments are easily modified to deal with the case in which some of the 
data samples are missing. In one embodiment, the calculations are performed only over 
those blocks where there are no missing data. In an alternative embodiment, missing data 
are set to zero, and the calculations are performed only over those blocks in which there 
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are nonzero data. The complexity calculations in process at block 200 of FIG. IB can be 
suitably modified to take this variation into account. 
A TOP-DOWN APPROACH 

Certain embodiments use a bottom-up approach in which all block-level correlations are 

5 computed first by the process at block 400 of FIG. IB, and are then used to calculate 
the approximate correlation sums D{f,a) for all candidate / and a. In an alternative 
embodiment, one may start by considering a limited set of candidate / and a, and calculate 
D{f,a) only for the chosen candidates. This limited set may consist of just one or multiple 
candidates. The calculation of D{f,a) for some candidate {f,a) requires certain quantities 

LO ui,h{cr) to be made available. These quantities can be calculated as needed. Furthermore, 
once a quantity ui,h{cr) is calculated, it is stored in memory, so that if it is ever needed in the 
future it need not be recomputed. This variation is useful in a situation where certain pairs 
{f,a) are examined first. If the ambiguity value associated with some pair (/,cr) turns out 
to be above some threshold, the process can stop. Otherwise, it can continue to examine a 

15 larger set of pairs. 

A SEQUENTIAL APPROACH 

Certain embodiments commit to processing the full set of data. In an alternative em- 
bodiment, a smaller set of data is examined first. This can be accomplished, for example, 
by keeping the same block length A^o but using a smaller number L' of blocks, so that only 

20 L'No<N of the data are utilized. If the ambiguity value associated with some pair (/, a ) is 
found to exceed a certain threshold, the algorithm is deemed to have successfully acquired 
the signal and the computation can stop. Otherwise, a larger number of blocks is fed to the 
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algorithm, and the algorithm repeats on the expanded set of blocks and data. Any inter- 
mediate results that have been calculated when using a smaller number of blocks are saved 
in memory so that they can be reused later. When dealing with a smaller set of data^ the 
grid-spacing of the frequencies to be searched can be taken to be coarser, with an attendant 
5 reduction of the computational effort. If the number of blocks is subsequently increased, 
the frequency grid is refined. With such an approach, the computational complexity of the 
overall method will be different, and suitable changes have to be made in the process used 
to select the various algorithm parameters. 
A HIERARCHICAL VERSION 
ilO In certain embodiments, there are only two frequency grids and the data are broken into 

\ a single set of blocks. In alternative embodiments, the same method is used in a hierarchical 
manner. In one embodiment, the blocks are broken down into subblocks and the calculation 
of the various block-level correlations is approximated by calculating correlation sums for 
each subblock and then merging them together to synthesize block-level correlations. This 

15 decomposition can be continued recursively by introducing an arbitrary number of levels. 

In the more general embodiment, a hierarchy of levels is used. With each level r there is 
an associated decomposition of the data into blocks, and an associated grid :7> of frequencies. 
For each level r, the grid of frequencies need not be regularly spaced. The only constraint 
is on the block sizes at the different levels. A block at a higher level r must be the union of 

20 smaller blocks at lower levels < r. One then proceeds as follows. Starting from the lowest 
level, call it level 1, calculate correlation sums a) for the various blocks i and for the 
various frequencies / in the set of frequencies J^i associated with that level. For a higher 
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level r > 1, and for all the blocks and all the frequencies in the set associated with that 
level, calculate approximate correlation sums, as a linear combination of the correlation sums 
calculated from a lower level. In particular, the calculation of a correlation sum Di{fya) at 
a certain frequency f £ J^r and for a certain r-level block makes use of the correlation 
sums Dj{f\ a) for those (r — l)-level blocks j that comprise block i, and for those frequencies 
/' e jF^_i that are closest to /. 
EXPLORING A SMALL FREQUENCY RANGE 

In some situations, one may want to evaluate approximate correlation sums a) for a 
small range of frequencies, or even at a single frequency. This is accomplished by forming a 
suitable linear combination of block-level correlations that are already available at frequencies 
in the vicinity of the frequency / of interest. An embodiment of this kind is of interest in 
various contexts. One particular context is when the signal has already been acquired, a 
moderately accurate estimate / of the true carrier frequency is available, and a search for a 
more accurate estimate by exploring the vicinity of / is desired. 

HARDWARE OVERVIEW 
An embodiment of the invention may be implemented using a computer system that 
includes a processor for processing information. The Computer system also includes a main 
memory, such as a random access memory (RAM) or other dynamic storage device, coupled 
to a bus for storing information and instructions to be executed by the processor. The main 
memory also may be used for storing temporary variables or other intermediate information 
during execution of instructions to be executed by the processor. The computer system 
further includes a read only memory (ROM) or other static storage device coupled to the 
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bus for storing static information and instructions for the processor. A storage device, such 
as a magnetic disk or optical disk, is provided and coupled to the bus for storing infornaation 
and instructions. 

The invention is related to the use of the computer system for implementing the tech- 
5 niques described herein. According to one embodiment of the invention, those techniques 
are implemented by the computer system in response to the processor executing one or more 
sequences of one or more instructions contained in main memory. Such instructions may be 
read into the main memory from another computer-readable medium, such as the storage 

tfl device. Execution of the sequences of instructions contained in the main memory causes 

ffi 

the processor to perform the process steps described herein. One or more processors in a 

%%$ 

m multi-processing arrangement may also be employed to execute the sequences of instructions 
- contained in the main memory. In alternative embodiments, hard-wired circuitry may be 
^: used in place of or in combination with software instructions to implement the invention. 
Thus, embodiments of the invention are not limited to any specific combination of hardware 
15 circuitry and software. 

The term "computer-readable medium" as used herein refers to any medium that par- 
ticipates in providing instructions to the processor for execution. Such a medium may take 
many forms, including but not limited to, non- volatile media, volatile media, and transmis- 
sion media. Non-volatile media includes, for example, optical or magnetic disks, such as 
20 the storage device. Volatile media includes dynamic memory, such as the main memory. 
Transmission media includes coaxial cables, copper wire and fiber optics, including the wires 
that comprise the bus. Transmission media can also take the form of acoustic or light waves, 
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such as those generated during radio wave and infrared data communications. 

Common forms of computer-readable media include, for example, a floppy disk, a flexible 
disk, hard disk, magnetic tape, or any other magnetic medium, a CD-ROM, any other optical 
medium, punch cards, paper tape, any other physical medium with patterns of holes, a RAM, 
5 a PROM, and EPROM, a FLASH-EPROM, any other memory chip or cartridge, a carrier 
wave as described hereinafter, or any other medium from which a computer can read. 

Various forms of computer readable media may be involved in carrying one or more 
sequences of one or more instructions to the processor for execution. For example, the 
instructions may initially be carried on a magnetic disk of a remote computer. The remote 
iXO computer can load the instructions into its dynamic memory and send the instructions over 
I a telephone line using a modem. A modem local to computer system can receive the data on 
the telephone line and use an infrared transmitter to convert the data to an infrared signal. 
An infrared detector coupled to the bus can receive the data carried in the infrared signal 
and place the data on the bus. The bus carries the data to the main memory, from which 
15 the processor retrieves and executes the instructions. The instructions received by the main 
memory may optionally be stored on the storage device either before or after execution by 
the processor. 

The computer system also includes a communication interface coupled to the bus. The 
communication interface provides a two-way data communication coupling to a network link 
20 that is connected to a local network. For example, the communication interface may be an 
integrated services digital network (ISDN) card or a modem to provide a data communication 
connection to a corresponding type of telephone line. As another example, the communi- 
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cation interface may be a local area network (LAN) card to provide a data communication 
connection to a compatible LAN. Wireless links may also be implemented. In any such 
implementation, the communication interface sends and receives electrical, electromagnetic 
or optical signals that carry digital data streams representing various types of information. 

The network link typically provides data communication through one or more networks to 
other data devices. For example, the network link may provide a connection through the local 
network to a host computer or to data equipment operated by an Internet Service Provider 
(ISP). The ISP in turn provides data communication services through the worldwide packet 
data communication network now commonly referred to as the "Internet" . The local network 
and the Internet both use electrical, electromagnetic or optical signals that carry digital data 
streams. The signals through the various networks and the signals on the network link and 
through the communication interface, which carry the digital data to and from the computer 
system, are exemplary forms of carrier waves transporting the information. 

The computer system can send messages and receive data, including program code, 
through the network(s), the network link and the communication interface. In the Inter- 
net example, a server might transmit a requested code for an application program through 
the Internet, the ISP, the local network and the communication interface. In accordance 
with the invention, one such downloaded application implements the techniques described 
herein. 

The received code may be executed by the processor as it is received, and/or stored in 
the storage device, or other non-volatile storage for later execution. In this manner, the 
computer system may obtain application code in the form of a carrier wave. 
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In the foregoing specification, the invention has been described with reference to specific 
embodiments thereof. It will, however, be evident that various modifications and changes 
may be made thereto without departing from the broader spirit and scope of the invention. 
The specification and drawings are, accordingly, to be regarded in an illustrative rather than 
a restrictive sense. 
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